Cooling dynamics of a dilute gas of inelastic rods: a many particle simulation 
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We present results of simulations for a dilute gas of inelastically colliding particles. Collisions are 
modelled as a stochastic process, which on average decreases the translational energy (cooling), 
but allows for fluctuations in the transfer of energy to internal vibrations. We show that these 
fluctuations are strong enough to suppress inelastic collapse. This allows us to study large systems 
for long times in the truely inelastic regime. During the cooling stage we observe complex cluster 
dynamics, as large clusters of particles form, collide and merge or dissolve. Typical clusters are 
found to survive long enough to establish local equilibrium within a cluster, but not among different 
clusters. We extend the model to include net dissipation of energy by damping of the internal 
vibrations. Inelatic collapse is avoided also in this case but in contrast to the conservative system the 
translational energy decays according to the mean field scaling law, E{t) oc t~ 2 , for asymptotically 
long times. 
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I. INTRODUCTION 

In a recent paper Q, hereafter referred to as I, we dis- 
cussed the properties of inelastic two particle collisions, 
starting from a Hamiltonian model for one dimensional 
elastic rods. Within this model, the coefficient of resti- 
tution e does not only depend on the relative velocity of 
the colliding particles, but in addition becomes a stochas- 
tic quantity, depending on the state of excitation of the 
internal vibrations. Here we extend the analysis to a dis- 
cussion of the many body dynamics of a one-dimensional 
gas of granular particles, modelled as elastic rods. We 
concentrate on dilute granular systems in the "grain in- 
ertia" regime, where two particle collisions dominate the 
dynamics. It was shown in I that successive collisions are 
to a very good approximation uncorrelated, so that the 
many body dynamics is a random Markov process. Con- 
sequently, the collisions are simulated by a Monte Carlo 
algorithm: velocities are updated with a random coeffi- 
cient of restitution, drawn from the appropriate proba- 
bility distribution. In between collision events, particles 
move freely like in an event driven algorithm. 

We focus here on the cooling properties of a large sys- 
tem (10000 particles) in the inelastic regime and refer to 
cooling as the decay of translational energy with time. 
We observe the evolution of spatial structures, without 
running into problems with inelastic collapse, which is 
always avoided by the algorithm. The most prominent 
spatial structures are large clusters of particles, which are 
seen to form and decay by colliding with other clusters. 
The velocity distribution within a cluster is to a good ap- 
proximation Maxwellian, whereas the global velocity dis- 
tribution shows significant deviations from Maxwellian, 
indicating that local equilibrium has been established 
within a cluster, but not among different clusters. 

The model is extended to include net dissipation, i.e. 
irreversible energy loss, in a phenomenological way by 
simply introducing a single relaxation time for the decay 



of the energy of internal vibrations. The final state of this 
model is one big cluster with all particles at rest. The 
dynamics with dissipation resembles a deterministic sys- 
tem (i. e. with constant coefficient of restitution) as long 
as no inelastic collapse is threatening. When the col- 
lision frequency increases dramatically, then ultimately 
the time between collisions will become smaller than the 
decay time for the internal energy, so that the vibrations 
no longer decay in between collisions. Then the model 
effectively reduces to the above stochastic one without 
dissipation and no inelastic collapse occurs. The kinetic 
energy of translation follows on average the mean field 
scaling law E(t) oc t~ 2 . 

Several groups have simulated onedimensional granu- 
lar media using event driven algorithms with constant 
coefficient of restitution. For (1 — e) > const. /N a di- 
vergence of the collision frequency in finite time, i.e. in- 
elastic collapse is observed. This leads to a breakdown of 
the algorithm and one either has to restrict oneself to the 
quasielastic regime, where e is sufficiently close to one or 
additional assumptions about the dynamics of clusters 
have to be made. Bernu and Mazighi || investigate a 
column of beads colliding with a wall. McNamara and 
Young jj] and Sela and Goldhirsch ||] discuss the cooling 
dynamics of a granular gas in the quasielastic regime. 
They observe the evolution of spatial structures and a 
bimodal velocity distribution. The critical wavelength of 
the instability is related to the minimum number of par- 
ticles for inelstic collapse to occur, given a fixed value of 
e. Clement et al. |J and Luding et al. study a vertical 
column of beads in a gravitational field with a vibrating 
bottom plate. For e close to one they observe a fluidiza- 
tion transition, whereas for t « 1 a bifurcation scenario 
is seen to take place. The latter has also been observed 
by Luck et al. (g) for a single bead on a vibrating plane. 

All of the above simulations use a coefficient of restitu- 
tion which is independent of the impact velocity whereas 
experiments on ice spheres reveal a velocity dependence 
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of e [0|. There have also been several attempts to calcu- 
late the velocity dependence of the coefficient of restitu- 
tion by extending the static theory of Hertz |[o| to vis- 
coelastic behaviour. One either assumes a phenomeno- 
logical damping term jnj in the equation of motion for 
the deformation or alternatively uses a quasistatic ap- 
proximation for low relative impact velocities. As 
a result of either approximation the coefficient of resti- 
tution becomes velocity dependent. Simulations of large 
systems with strongly inelastic collisions have been per- 
formed with this model |Q. 

Another approach is based on phcnomcnological wave 
theory. Here one assumes that two colliding bodies do 
not vibrate before collision and that the impact triggers 
a travelling clastic wave in both of them. For one dimen- 
sional rods this ansatz yields [ )l5| e = h/h, independent 
of the relative velocity of the colliding particles. Here l\ 
(h) denotes the length of the shorter (longer) rod. As 
shown in I, these results are contained in our model. 

Our paper is organized as follows. In the next sec- 
tion we review the model of I, discuss the probability 
distribution of e and define the algorithm for the many 
body dynamics. Results of simulations are presented in 
We first discuss global quantities, like the time 
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sec. 

decay of the kinetic energy and the total number of colli- 
sions as a function of time. Subsequently we analyse the 
local structure with help of the paircorrelation function 
and discuss the formation and decay of particle clusters, 
as well as the distribution of the particles' velocities. In 
sec. 



IV the model with dissipation of vibrational energy 
is introduced. Finally in sec. ^ we summarize our results 
and give an outlook to forthcoming work. 



II. MARKOVIAN DYNAMICS OF INELASTIC 
RODS 

We first review the Hamiltonian model of I and sum- 
marize the properties of two particle collisions. We then 
show that the transition probabilities of the resulting 
Markov process obey detailed balance and introduce the 
algorithm for the dynamics of the many-body system. 



A. Two-particle collisions 

Our starting point are the Hamiltonian equations of 
motion of a system of N elastic rods of homogeneous mass 
density. The particles are placed on a ring of circumfer- 
ence L. Each rod is characterized by its length total 
mass rrii and centre of mass position Ri(t). Its vibra- 
tional excitations are described by N moc i normal modes 
q\(y = l,...,N mo d) of wavenumber fcj ; „ = nu/li and 
frequency u)i >v = cki M . The only important material 
parameter for our model is the sound velocity c. We 
model collisions of the rods by a short range repulsive 
potential V(r) = _Bcxp(— ar), which depends on the 



momentary end-to-end distance r between the colliding 
rods, thus coupling translational and vibrational degrees 
of freedom. We shall be interested in the hard core limit, 
which can be achieved by letting a — > 00 The total 
Hamiltonian of our model reads: 



n 



Hbath {pj" \ qf ] } + Htr {Pi} + n int {i?,, 9l M } 



N N mod [ (»2 
ST ST ) P i , 2 Q; 



e5-+ m 



2=1 U=l 

N-l 



2m i 



2m 



2^ Be V 



Here Ri+i,i = Ri+i — Ri — (h+i — U)/2 is the end-to-end 
distance of two undeformed neighbouring rods and Pi(t) 
and Pi(t) denote the conjugate momenta for the centre 
of mass and the amplitude of vibration. The first term, 
'Hbath, models the internal vibrations, the second term, 
Ti.tr, the translational motion, and the third term, Hint, 
the interaction between the rods. 

In I we have analysed the statistical properties of 
two-particle collisions by numerically integrating the full 
Hamiltonian dynamics for the case N = 2 with length 
ratio 7 = h/h- The main results are the following. 
Equipartition among the vibrational states of a rod is 
achieved fast (after about five collisions) as compared to 
the relaxation of the translational velocity, which hap- 
pens on a time scale of about 80 collisions. The coeffi- 
cient of restitution of two successive collisions is to a very 
good approximation uncorrelated. Based on these obser- 
vations a simplified description was achieved in I by inte- 
grating out the internal degrees of freedom, which can be 
done exactly for a single two particle collision. One is left 
with an effective equation of motion for the rescaled rel- 
ative velocity of the two rods w(t) = R.2,i(t)/ R2,i{0) — 1. 
Here r = ct/l and I = 2l\l2j(l\ + h) denotes an effective 
length which is always chosen much larger than the range 
of the potential, i.e. al ^ 1. Based on the observation 
of fast equipartition among the vibrational states, these 
are modelled by a thermalized bath, characterized by a 
temperature Tb = Ebath I ^ mod, where the vibrational 
energy of a rod is given by the sum of the energies of the 
individual modes: Ebath — Y^=\ d E v - Thus q"(0) and 
Pi(0) are taken as independent, canonically distributed 
random variables 
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"The constant B is arbitrary, it can be absorbed by rescaling 
time t —> t\fi3 and frequencies uj — > uj/s/B. 
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Under these assumptions the relative velocity w(t) obeys 
a stochastic equation of motion 

—w(t) = — cxp { k(t — To — w(t) — 
or k i 

2 Co n 

£5>(T-nr,) + g(T)) , (3) 

i=l n=l ' 

where k = — s ^Rn,i(0) and r = ai? 2 ,i(0)/K. The coeffi- 
cient of restitution is given by 



e = lim w(t) — 1 



(4) 



and thus becomes a stochastic variable, depending on the 
state of the vibrational bath before the collision. q(r) is 
a Gaussian random noise with zero mean and covariance 



C q (r) = (q(r')q(r'+r)) 




(5) 



n=l 



Here, 9 denotes the Heaviside step function and E tr = 
A i (^2,i(0)) 2 /2 is the translational energy of the colliding 
rods in their center of mass frame of reference. The Ti 
which appear in eqs. (||) and (|^) are determined by the 
length ratio 7 of the rods according to T± = 1 + 7 and 
r 2 = 1 + ^ and /i = TO 1 r7i 2 /(m 1 + m 2 ) is the effective 
mass. The stochastic process {q(r)} is simply related to 
two periodic brownian bridge processes with periods 
Ti and r 2 respectively 

In the hard core limit (k — > 00) the stochastic equation 
for the rescaled relative velocity w(t) (eq. (||)) can be 
solved by saddle point methods, yielding 

w(t) = max (0, /(r)) , where 

f( T ) = o^xj T ' - T o - W ( T ' ~ ^ r ') + ? ( T ') f ■ ( 6 ) 

The duration of the collision as well as the final velocity 
are stochastic variables. The collision is ended when the 
memory terms in eq. ([}]) overcompensate the gain from 
the other terms t' — tq + q(r'), which are on average 
increasing. 



B. Transition Probability 



collisions. During a collision E tr changes to a new value 
E' tr — E tr e 2 . The probability for this transition is de- 
termined by the probability density for the coefficient of 
restitution Pfj{c) according to 



t (E tr -> E' tr ) = 
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(Tg here denotes the bath temperature of both rods, un- 
der the assumption that the temperatures are equal. If 
that is not the case, one would have to replace the index 
f3 by two indices [3\ and Here, we use only one index 
for notational simplicity.) 

Changes in the bath temperature are not independent, 
but determined by energy conservation: 
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The stationary state of the Markov process is known: 
after cooling, the system of two particles, each equipped 
with an internal bath, evolves into a stationary state with 
a Boltzmann distribution for Et r 
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where E to t is 



with the bath temperature 
the total energy of the system. 

It can be proven that this collision process obeys 
detailed balance, which gives the following relation for 
p 13(e) (this relation only holds if the temperature of both 
rods is indeed equal, i. e. /3i = /3 2 = [3): 



p f3 (e)e-f s =p e zJ-)e-^ 



(10) 



When the temperatures of the baths of oscillators are 
zercjj], i. e. q{r) = for all r, the collision is deterministic. 
For this case, the coefficient of restitution is equal to 7, 
the ratio of the lengths of the rods (see I and |15|). Thus 
in the limit of small temperatures, pp(£) should approach 
a 5- Function centered around 7. 

At large temperatures, on the other hand, simulations 
suggest the (quite sensible) result that pp(e) is a uniform 
distribution, i. e. all possible e's are equally probable. 

It can be shown from equation (^|) that for 7 = 1 the 
collision is always deterministic, i. e. e = 1 for any re- 
alisation of the stochastic process q(r). This interesting 
result will have implications on our setup of the simula- 
tions (see sec. II C). 



The results of the preceding section are interpreted as 
a Markov process in discrete time, which accounts for 
transitions of the translational energy upon successive 



'Actually, it is sufficient that only the temperature of the 
longer rod is zero. 
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C. Algorithm 



We now consider the dynamic evolution of N particles 



on a ring of circumference L + X^iLi h- L is thus the 
total length of the interparticle spacings. For the fol- 
lowing arguments the actual lengths of the particles are 
unimportant because the point in time when a collision 
occurs depends only on the end-to-end distance and the 
outcome of a collision depends only on the length ratio. 
In order to keep the notation simple we map the system 
to an equivalent one consisting of N point particles on 
a ring of circumference L. Each particle is characterized 
by its position Ri(t), its velocity Ri(t), the temperature 

of its internal bath Tg\t). The N mo d internal modes of 
one rod are represented by one degree of freedom only, 

namely Tg' = Y^v=i d l^mod- The rods are assigned 
alternating lengths such that the ratio 7 for each collision 
has a fixed value, in our case 0.8. The ratio of masses is 
also given by 7, assuming the same homogeneous mass 
density for both kinds of rods. We choose rods of al- 
ternating le ngth because due to the result given at the 
end of sec. II B a length ratio of 7 = 1 implies e = 1 
always which would correspond simply to a standard one 
dimensional hard sphere gas. 

The model we use is a hybrid of an event driven al- 
gorithm and a Monte Carlo simulation. The particles 
move freely in between collisions, as in event driven al- 
gorithms. When two particles collide their states are up- 
dated stochastically, according to the distribution of e. 

It is convenient to introduce dimensionless variables 
X{ = RiN/L and v% = Ri\/ h/Tq. Tq serves as an energy 
scale and will be identified with the homogeneous initial 
granular temperature of the many particle system. Time 
is measured in units of L^J fi/T /N. 

For the algorithm we only need relative distances and 
velocities 
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The algorithm is defined by iteration of the following 
steps: 

1. Calculate the time difference At for the next colli- 
sion to take place: 



( AXi 

At = mm — - — 

{i|A-u;<0} V Avi 



(13) 



The pair of particles which is going to collide next 
is denoted by (io,io + 1). 

2. The relative distances of all particles are updated 
according to: 



For the designated pair (io,io + 1) we obtain 
Ax l0 {t + At) = 0. 

3. The kinetic energy of relative motion of the pair 
(io, io + 1) as well as the mean local bath tempera- 
ture are calculated according to E tr — Avf g /2 and 

T B = (T£ o) +T^ 0+1) )/2. Subsequently, a random 
value of e is chosen from the probability distribution 
Pp(e), presently calculated by numerically solving 
eq. (g) and applying eq. (g). 

4. The bath temperatures and relative velocities are 
updated 



T ( z°> (t + At) = Tg 0+1 > (t + At) 
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(15) 



(17) 



Axi(t + At) = Axi(t) + Avi(t)At. 



(14) 



5. Continue with step 1. 



III. SIMULATIONS 

Many-body simulations using the above algorithm 
have been performed to study the cooling dynamics of 
the system. More precisely, we focus here on the inter- 
mediate range of timescales where equipartition among 
the internal modes has already been achieved and the fi- 
nal equilibrium state of equipartition among all degrees 
of freedom is not yet reached. As shown below this time 
range extends over several orders of magnitude. 

We assume that two particle collisions dominate the 
dynamic evolution of the system. This is justified for a 
dilute granular gas. The typical time of interaction in 
our model is given by ti nt = 2//c, i.e. the time a sig- 
nal needs to travel back and forth on a rod. Hence in 
principal, two colliding rods can interact with a third 
one. This will be highly unlikely, as long as the time be- 
tween collisions is much longer than U nt . This requires 
2l/c < L/(N\R i+ i <i \). So either the length of the rods 
has to be chosen sufficiently small as compared to the 
mean distance L/N or the initial velocities should 
be small compared to the velocity of sound. The latter 
is a material parameter and can have quite high values 
for hard materials (e.g. for steel, c ~ 10 4 m/s), favour- 
ing short interaction times. In a standard event-driven 
simulation inelastic collapse occurs when the number of 
particles is sufficiently large, resulting in a diverging col- 
lision frequency. This would clearly violate the condition 
that the time between two collisions is long compared 
to Unt- However, since our algorithm avoids the inelastic 
collapse, as will be discussed below, we will still make use 
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of the assumption that three or more particle collisions 
will not be important. 

A system of 10000 particles has been simulated and 
for most part a length ratio 7 = 0.8 has been used. 
We start from a spatially homogeneous distribution of 
particles having a maxwellian velocity distribution with 
(?f-u|) = 5 (i = 1(2) stands for the shorter (longer) 
species of rods). We use N mod = 1000 vibrational modes 

per particle. Initially the internal bath temperature T$ 
of each particle is set to 0. 

Our simulations were performed on a cluster of Linux 
workstations with Pentium processors. The longest runs 
took about three weeks of computer time. 



A. Global quantities 

Kinetic energy 

The time development of the total kinetic energy 
which is given by 



N Ti 

E k in = ~J V * 
i=l 



(19) 



in our rescaled units, is shown in fig. [I], in comparison 
with results for the deterministic model with constant 
e. For small times the curves for the deterministic and 
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FIG. 1. Reduction of the kinetic energy per particle as a 
function of time. The curve for the deterministic coefficient 
of restitution breaks off because an inelastic collapse occured. 
The dot-dashed line shows the average energy per vibrational 
mode for the 50 particle run and illustrates that equipartition 
holds in the stationary state. 

the stochastic dynamics are rather similar. In the initial 
stage very little energy is stored in the internal modes 
and hence the coefficient of restitution is approximately 
given by the deterministic value. However the determin- 
istic dynamics runs very quickly into the inelastic col- 
lapse, as can be seen from the total number of collisions 



which is shown as a function of time in fig. When this 
happens, the simulation gets stuck so that the curve for 
the kinetic energy breaks off in fig. [l]. The stochastic dy- 
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FIG. 2. Number of coilisions as a function of time. The 
inset shows a comparison of the deterministic model (dashed 
line) and the stochastic model (solid line). (The units on the 
axes of the inset are the same as on the regular axes.) The 
deterministic model quickly runs into the inelastic collapse, 
seen by the diverging number of collisions. The dotted lines 
shows the theoretical number of collisions as a function of time 
in the stationary state according to eq. ( ^l| ) for the 10000 and 
the 50 particle runs. The data for the 50 particie run has been 
scaled by a factor of 100 in order to fit on the graph. 

namics shows completely different behaviour: The kinetic 
energy continues to decrease until equilibrium is reached, 
where Ekin continues to fluctuate around the stationary 
value which is given by Ej£% = E kin (t = 0)/(2N mod + l). 
(The final state has not quite been reached for the 10000 
particle run in the time interval that is shown in fig. [j].) 

The final state of our stochastic model is a consequence 
of the idealised assumption that the total system be con- 
servative. In a more realistic model of granular media one 
expects the particles to be at rest in the final state. In 
sec. IV we shall present a phenomenological extension of 



our model, which does take into account energy dissipa- 
tion of the microscopic degrees of freedom. In the final 
state of the model with energy dissipation all particles 
are at rest. 

Collision rate 

Simple mean field arguments || have been used to de- 
rive scaling laws for the time evolution of kinetic energy 
and collision rate. One assumes that the particle veloc- 
ities are uncorrelated and Gaussian distributed. For a 
constant coefficient of restitution one obtains Ekin (t) ~ 
t~ 2 and N co u ~ Int. Neither scaling law fits our data, as 
can be seen from figs, [j] and |[ McNamara and Young 
H have already pointed out that the mean field scaling 
laws are only applicable in the quasi elastic regime, where 
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no inelastic collapse occurs. Otherwise the assumption 
of uncorrected Gaussian velocities breaks down. In the 
stochastic model we have additional fluctuations of the 
coefficient of restitution, which invalidate the derivation 
of the above scaling laws. Hence it is no surprise that the 
data disagree with these relations. 

The rate of collisions becomes constant as the station- 
ary state is approached, as can be seen from fig.^L The 
average collision rate is given by N C oii = NAv/ (2 Ax). In 
the stationary state the velocities are indeed uncorrelated 
Gaussian variables, distributed according to 
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(20) 



where j — 1(2) again stands for the shorter (longer) type 
of rods. We assume Ax — 1 and perform the average 
over velocities to obtain 



i\jstat 
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2y/TT(2N mod + 1) 



(21) 



This result is in very good agreement with the simula- 
tions of the 50 particle system in the stationary state 
(see fig. H). The 10000 particle system is also approach- 
ing the correct value as it gets closer to the stationary 
state. 



B. Local quantities 



Particle density 
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FIG. 3. Time evolution of the particle density. Dark re- 
gions indicate high density. 



8000 



6000 - 



4000 - 



2000 





_ 



10' 



10 



10" 



FIG. 4. The same as fig. |ona logarithmic time scale 



As is well known, inelastic particles without internal 
structure tend to cluster in one dimension. This clus- 
tering leads to a break down (inelastic collapse) of the 
system if e is constant and less than a critical value that 
depends solely on the number of particles. In our model 
particles have internal degrees of freedom and the trans- 
lational energy is not completely lost in a collision but is 
stored in the internal vibrations and can be transferred 
back to the translation al m otion. Due to the proper- 
ties of pp(e) (see section |ll B|) , the probability for this to 
happen gets larger as the translational energy decreases. 
Therefore clusters do form but dissolve after a while and 
no inelastic collapse takes place. 

We start again from a spatially homogeneous distribu- 
tion of particles and analyse evolving spatial structures 
with help of a coarse grained density p. We divide the 
total length of the ring into a hundred bins and count 
the number of particles in each bin. The coarse grained 
density is defined as the actual number of particles in 
each bin divided by the average. 

The time evolution of p is shown in fig. |^ on a linear 
time scale and in fig. |ona logarithmic time scale. Sev- 
eral phases in the cooling process can be identified. First, 
the particles start to form clusters and voids as they lose 
kinetic energy in collisions (initially, when Tg is small 
compared to the translational energy, the coefficient of 
restitution is always close to 7). After these clusters have 
formed, one observes collisions of clusters, forming larger 
clusters. Simultaneously the dissolution of clusters starts 
to set in, the remains being sent outwards to join neigh- 
bouring clusters. The biggest clusters and voids are seen 
to survive for times of order 10 4 . This complex interac- 
tion of forming and dissolving clusters continues with a 
clear tendency to form fewer and larger clusters. Finally 
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these large clusters dissolve to establish the equilibrium 
state, i. e. equipartition among all degrees of freedom. 

For 10000 particles it takes a time of order 10 5 un- 
til the cooling dynamics is finished and the equilibrium 
state for the larger system is reached, whereas for 50 
particles it takes only a time of order 10 3 (see fig. [l]). 
The equilibrium state is reached only after the formation 
and dissolution of essentially one final large cluster. In a 
smaller system the end of the cascade of clusters of in- 
creasing size is reached earlier, simply because there are 
fewer particles. 

Phase space 

The complete information about the state of the sys- 
tem at time t is contained in a phase space plot, as shown 
in fig. H Within a cluster of particles we expect frequent 
collisions and hence an effective transfer of kinetic energy 
to internal vibrations. Frequently regions of high average 
density are characterized by particle velocities centered 
around zero. However, we also observe clusters with an 
average nonzero velocity, resulting at a later time in col- 
lisions of clusters. One such collision of two clusters can 
be traced in fig. || around x ~ 2500. In fig ||a (a snap- 
shot taken at t = 14000) one observes two clusters both 
with nonzero average velocity moving towards each other, 
wheras in fig ||b (taken at t = 20000) the clusters have 
collided and formed a larger one. 

We also see around x ~ 1000 the occurrence of a stripe- 
shaped fluctuation in the phase-space plot. This type 
of fluctuation has already been observed and discussed 
by McNamara and Young || and Sela and Goldhirsch 
||. It gives rise to the formation of clusters out of an 
initially homogeneous region. Thus figure ^| shows that 
the dynamics of our system are indeed rather complex 
as formation, movement, interaction and dissolution of 
clusters all happen simultaneously. 

Local kinetic energy 

It is interesting to see how the kinetic energy is spa- 
tially distributed. We define a coarse grained kinetic en- 
ergy density similar to the coarse grained density by sum- 
ming the kinetic energies of all particles inside a bin and 
dividing by the number of particles in the bin. 

One might be tempted to conjecture that the local ki- 
netic energy is in some way correlated to the clustering 
because most collisions occur within the clusters. Fig. g 
(as an example) reveals, however, that this is generally 
not the case: although the kinetic energy shows some 
structure there is no visible correlation to the density, 
not even in a state as the one shown in fig. Iq, where all 
the particles are extremely clustered. Fig. g is a snap- 
shot of the system at time t — 100343 (cf. figs. || and||). 

Velocity distribution 

In the cooling stage, the system is still far from equilib- 
rium, so that the velocity distribution of the particles is 
not expected to be a Maxwell distribution. It is therefore 
interesting to test what kind of distribution the velocities 
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FIG. 5. Phase space plot of the system at two different 
times. 



really follow. 

Data analysis shows that the velocity distribution of 
all particles is indeed not a Gaussian distribution (see 
fig. 0). There are relatively large deviations especially 
near the maximum of the curve. If one restricts the data 
analysis to only those particles inside a single cluster, 
however, one finds that the velocity distribution of these 
is to a much better degree gaussian, considering that 
there are only about l/10th of the total number of parti- 
cles in the cluster. This can be well understood because 
there are many collisions between particles inside a clus- 
ter and thus a local equilibrium is reached, resulting in a 
Maxwellian distribution. On the other hand the velocity 
distribution of all particles reflects the velocity distribu- 
tion of the clusters. As long as the complicated process 
of forming and dissolution of clusters is underway, the 
clusters are naturally far away from equilibrium. This 
leads to the observed deviations from the gaussian curve. 

Since our system is far away from the quasielastic limit, 
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FIG. 6. Comparison of the local kinetic energy (dashed 
line) and particle density (solid line). 



we see quite a different velocity distribution than Mac- 
Namara and Young [Q , who simulated a one dimensional 
system of quasielastic particles. They observed a bimodal 
velocity distribution because the particles tend to concen- 
trate on the upper and lower edges of a band in a phase 
space plot similar to fig. |^. In our simulation, the situ- 
ation is much more complex because of the formation of 
many clusters, each with its own velocity distribution. 

Correlation function 

The inelasticity of collisions leads to a clustering of 
particles, as can be seen in figs. | and H Williams |J] 
has described a one dimensional system of individually 
heated granular particles. He found that the pair correla- 
tion function, defined by g(x) = j^—^ 5{x— \xi — Xj\) 
of the system in the steady state approximately follows a 
power law. Here, we observe quite a different behaviour 
of the correlation function (see fig. ||) . Instead of showing 
a divergence at zero separation, it levels off to a plateau. 
The explanation for such a different behaviour lies in the 
mechanism of heating: When the particles are heated in- 
dividually, i. e. when they are driven by a random force, 
they will half of the time be kicked back in the direction 
of the particle that they last collided with. Thus there 
is some additional tendency for the particles to stick to- 
gether. In our model, however, the particles will only 
change their velocity when they collide, thus favouring 
larger distances. 

It should be noted that the correlation function in fig. 
|§| is not that of the steady state of our system but a snap- 
shot taken during the cooling process. The steady state 
of our model is trivial, implying a constant correlation 
function, g(x) = 1. 



FIG. 7. Velocity distribution of all particles (circles) and 
the particles inside one particular cluster (triangles) at time 
t — 14000. The cluster chosen for this curve is centered 
around x = 2000 (cf. figs. |, § and |a) 

IV. DAMPED INTERNAL MODES 

So far we have considered a conservative system, i. e. 
the total energy of translational motion and internal vi- 
brations is conserved. Such a model gives rise to a sta- 
tionary equilibrium state in which equipartition among 
all degrees of freedom holds, so that the translational 
momenta are of order 0(l/y/N mo d). To model granular 
media, one should take into account additional dissipa- 
tive mechanisms, which result in a decrease of the total 
energy so that the particles are truly at rest in the sta- 
tionary state. One such mechanism is black body radia- 
tion. 

A simple way to model this effect is to let the bath 
temperature of each particle decay in time. Hence we 
sugg est the following modification of the algorithm of sec. 
IIC. Inbetween collisions the particles move freely and 



their bath temperature decreases according to a simple 
exponential decay 

Thit) = T'h{U) exp (-(* - U)v) for t > U (22) 

Here tj denotes the instant of time when the last collision 
of particle i took place. The same decay frequency v is 
used for all particles. The updating of relative veloci- 
ties and bath temperatures in a collision is unchanged as 



[IC 



compared to sec. 

We expect that the effect of such a dissipative mecha- 
nism will strongly depend on the frequency of collisions 
as compared to v. If collisions are rather infrequent, then 
the decay will be effective and the bath of the particles 
will cool down in between collisions. The resulting dy- 
namics should resemble the deterministic case and hence 
one should observe a strong increase in the collision fre- 
quency, because the system is developing towards inelas- 
tic collapse. When this happens, the collision frequency 
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FIG. 8. The pair correlation function g(x) of the system at 
t = 14000. 



FIG. 9. Number of collisions as a function of time for the 
dissipative system. 



becomes comparable or even smaller than the decay rate 
v. In that case the internal modes can no longer re- 
lax in between collisions. In the limit of very high col- 
lision frequencies the bath temperatures are effectively 
non- decaying, so that one recovers the algorithm of sec. 
II C without any dissipation. Hence we expect to see the 
system develop towards inelastic collapse with a strong 
increase in collision frequency, followed by a period of 
time where the collision frequency levels off. 

These expectations are confirmed by numerical simu- 
lations of once again 10000 particles with v = 0.01. In 
fig. ||we show the total number of collisions as a function 
of time. One clearly observes rather sharp steps followed 
by smoother regions, as explained above. In fig. ^ we 
show the decrease in total kinetic energy as a function of 
time. Steep regions of the collision frequency correspond 
to steep regions in the energy plot because the frequent 
collisions among clustered particles draw much energy 
out of the system. 

An important question is the following: What is the 
stationary state of the system with dissipation and how 
long does the system need to relax to the stationary 
state? The final state should be one big cluster, with 
all particles at rest. As explained above, the dynamics 
with dissipation resemble the dynamics of a determinis- 
tic system as long as no inelastic collapse is at hand. For 
this reason it can be expected that the kinetic energy on 
the average follows the mean field result E ~ t~ 2 , occa- 
sionally disrupted by the occurence of a cluster, which 
is however quickly dissolved. Simulations show that this 
behaviour can indeed be observed, see fig. [□]. Since the 
above mentioned scaling law never permits the energy to 
become exactly in a finite time, the stationary state 
(which has energy 0) can also never be reached in finite 
time. 
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FIG. 10. Decrease of the total kinetic energy as a function 
of time for the dissipative system. 



V. CONCLUSION 

We have presented the results of simulations performed 
on a recently developed model for a one dimensional gran- 
ular medium. The model allows for an algorithm which 
is a hybrid of Monte Carlo and Event Driven and which 
avoids the inelastic collapse. Thus, we have been able 
to perform long simulations on a large system far away 
from the quasielastic limit. The model in its simplest 
form conserves energy: translational energy can be trans- 
ferred to internal vibrational modes of the particles and 
vice versa. Starting from a state with no internal modes 
excited, there is a long cooling regime, extending over 
several orders of magnitude in time, before the stationary 
state, characterized by equipartition of energy, is reached. 
The decrease of kinetic energy during this cooling stage 
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FIG. 11. Decrease of the kinetic energy as a function of 
time for 50 dissipative particles. The dashed line shows the 
E ~ t~ 2 relation. 

shows considerable deviations from the mean field result 
Ekin ~ t~ 2 - We have also observed a complex process of 
cluster forming, movement, interaction and dissolution. 
Inside the clusters we find that the particles are close to 
local equilibrium, which is indicated by the fact that a 
Maxwellian velocity distribution with in general nonzero 
mean velocity holds for particles in the cluster. 

The model has been extended to include net dissipa- 
tion of energy by exponential damping of the internal 
modes. In this case, the algorithm still shows no inelas- 
tic collapse. On the average, the decrease of the kinetic 
energy follows the mean field result but with considerable 
fluctuations due to the complex cluster dynamics which 
are still a feature of the model with energy dissipation. 

Thus our model, which is based on a microscopic mech- 
anism for the loss of translational energy during colli- 
sions, is well suited as a starting point for simulation and 
theoretical description of one dimensional granular me- 
dia. Unlike many other models, it makes use of an exact 
treatment of the collision dynamics of the colliding rods 
and hence offers a possible intuitive way of understand- 
ing the precise manner in which translational energy is 
removed from a granular system. Future work will use 
this model to investigate the properties of driven granular 
assemblies. In our model a specific mechanism - transfer 
of translational energy to internal vibrations - has been 
analysed to develop a microscopic basis for an effective 
coefficient of restitution. One may wonder which of our 
results depend on the particular mechanism. To study 
this question, we are presently investigating distributions 
Pfj(e) which are only restricted by detailed balance and 
not derived from a microscopic model. One may also 
try to extend our analysis to higher dimensional objects 
like disks or spheres. In the simplest geometry these ob- 
jects are colliding in a one dimensional tube, so that no 
tangential forces like coulomb friction have to be consid- 



ered. Our microscopic model is easily generalized to this 
case and allows for investigating how effectively energy of 
translation is transferred to elastic vibrations [ p0[ . The 
frequently used quasistatic approximation of Hertz im- 
plies no energy transfer at all. 
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